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The dynamical density-density correlation function is cal- 
culated for the one-dimensional, half-filled Hubbard model 
extended with nearest neighbor repulsion using the Lanczos- 
algorithm for finite size systems and analytically for large on 
site repulsion compared to hopping amplitudes. At the zone 
boundary an excitonic feature exists for any finite nearest 
neighbor repulsion and exhausts most of the spectral weight, 
even for parameters where no exciton is visible at zero mo- 
mentum. 



Materials such as SrCuO^ and Sr 2 Cu03 are believed 
to be quite well described by the one-dimensional half- 
filled Hubbard modetl. With the current level of experi- 
mental technology it is possible to measure dielectric re- 
sponse (which is directly related to the density response 
function) withpgood resolution in both momentum and 
energy transferal, which makes it of great interest to un- 
derstand the response expected from this simplest model 
of a Mott-Hubbard insulator. In this paper we present 
the combined results of analytic and numerical calcula- 
tions which clarify the nature of the dynamical density 
response function of the Hubbard model including also 
nearest-neighbor repulsion, which is a first step to the 
inclusion of the long-range Coulomb interaction. The 
results presented may serve as a indication of the ap- 
propriateness of the model to particular materials once 
experimental data become available. Previous work of 
Mori, Fukuyama and Imadau in this direction was based 
on an effective low-energy model whose relation to our 
extended Hubbard model is non-trivial. 

The extended Hubbard model is defined as 



Recr(w) = lim = Af(k,uj). 

y ' fc-0 4sin 2 (fc/2) V ; 



(3) 



n 



h.c. +U + V, 



(1) 



where hi = Tij^rij^ is the on-site Coulomb repul- 

sion, V = VJ2j n j n j+i i s the nearest-neighbor repulsion 
and rii = riii + nil IS the density at site i. 

The imaginary part of the density-density correlation 
function is given by 

Af(k,iu) = |(/K|GS)| 2 (5( W - E f + Eqs), (2) 
/ 

where |GS) is the half-filled ground state and |/) denotes 
a final state with energy Ef and momentum k = Pj. 
Note that the optical conductivity is directly related to 
the density-density correlation function through the con- 
tinuity equation (see e.g. Ref. [p|) 



For the case of small charge gap compared to bandwidth 
a(uj) has been calculated by Giamarchi and MillisEl us- 
ing bosonizatioa, and the large U case was studied by 
Gebhardt et alB 

We evaluate Eq. (||) by: a) direct diagonalization, 
which is limited to small clusters; b) for U ^S> t,V we per- 
form a canonical transformation of the type introduced 
by Harris and LangeH to get a a t J" like effective model. 
This model is then diagonalized for larger system sizes; 
c) finally, for the effective model in the thermodynamic 
limit an analytic expression is obtained-,using the factor- 
ized wave function of Ogata and Shibau. 

a. Exact diagonalization of Hubbard model: Using the 
standard Lanczos algorithm, it is straightforward to cal- 
culate Af(k,ui) for finite size clusters of L sites. In Fig. |l| 
we show a typical plot for relatively large U and inter- 
mediate V. There are several features to observe: the 
spectra consist of dominant peaks which are distributed 
over a frequency range from approximately U — V — At 
to U — V + 4t at small momentum transfer, narrowing to 
a single peak near U — V at the zone boundary. There 
are also less intense features which are barely visible in 
the figure whose weight increases slightly for smaller U . 
The number of large peaks at small momentum scales 
linearly with the system size, and is therefore most nat- 
urally interpreted as the finite-size precursor to a contin- 
uum absorption. Near the zone boundary however, for all 
accessible system sizes we see only a single peak when V 
is finite. We shall see that this is due to the appearance 
of an exciton below the continuum. A further feature of 
Af(k, ui) is the asymmetry in the weight distribution: the 
spectra are skewed toward lower frequency, which is more 
pronounced for either smaller U or greater V. Our main 
interest here is the region however we note that 
for V approaching U/2 the transition to the charge den- 
sity wave statety is directly seen as a drastic softening of 
the response at k = n. Apart from the features presented 
there is of course also the special case of k = where the 
response is simply a delta function at zero frequency. 

b. Effective model: In order to understand more 
clearly the above observations it is instructive to derive 
an effective model in the strong-coupling large U lim^ 
We closely follow the approach of Harris and Lange 
however, for details and notation we refer to Ref. [ [ft). 
We first define 



To = V 



1 



-t ^(1 - fti,s)cl )<7 c. +J)(7 (l - n i+S ,a) 

4,(5,(7 

T U = ~tYl ^^[A+sA 1 ~ (4) 



i.S.t 



and TLj/ = T^, where the subscript denotes the change in 
the eigenvalue of U induced, corresponding to the change 
in the number of doubly occupied sites, so that 



H ctf =U + f + -^ 



Tu,T- 



Here by O we denote the canonically transformed op- 
erators. The expression for TC c s at the operator level is 
quite long and complicated, here we rather give the effect 
of U~ 1 [Tjj, T-u] applied to the basis states: 
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(2t 2 /U)(\aa) 



(2tyu)(\T?) + \?T)) 

(t 2 /U)(\aar)-\aar)) 
(t 2 /U)(\raa)^\raa)) 

-(t 2 /U) (|oTt) + \<TfT)) 

-(t 2 /U) (\rfa) + \fra)) , 



(6) 



where a stands for spins and r for empty (e) or doubly 
occupied (d) states, and we use the convention where 
the creation operators for the states are ordered with 
increasing site index and in case of double occupancy we 
usec| T ct ; . 

In order to calculate the matrix elements describing the 
transition to the upper Hubbard band in Eq. (Q) , we must 
transform the electron density operator appropriatelyO, 
ni,u = n\if + nflj + . . ., where the leading term n^j = 
{T(j,hifi]/U reads 



S=±l 



(7) 



The next order correction will be considered later. 

At half-filling the ground state wave function of Ti e s 
is simply that of the Heisenberg model, and matrix ele- 
ments are with the states in the upper Hubbard band - 
i.e. with the states containing exactly one doubly occu- 



pied site. Then (/|»$|GS> = t/U(f\ E CT , 5 (c] 
)|GS) holds and Af(k,u>) simplifies to 



f c 



tf(k,u>) = £^4 sin 2 \ |</I(4A ct " 4<A CT )|GS) 



U 2 

xS{lo -E f + E GS )S kt 



Pt- 



(8) 



The density response for this effective model as deter- 
mined by exact diagonalization is shown in Fig. |l| as the 



dashed curve. We can see that the overall behavior fol- 
lows that of the Hubbard model, however there are sig- 
nificant deviations in the distribution of weights near the 
edges of the spectrum for the not extremely large U used. 
This agreement may be improved by including also the 
next to leading order correction to the density operator: 



~(2) 

n 3,V - Jj2 
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j+26,a C j+S,a C j+S,B- C j,a 



' 7/2 {pjcr^j+Sv C j+S<r C ja- 

(7,5 



./..t' , + ,v,t' j + .y.-i ' :--2S,a)] 

(9) 



where the final form is correct only when applied to a 
half-filled state with no double occupancy. The diago- 
nalization result including also these terms in the matrix 
elements is shown as the dotted line in Fig. [|. The agree- 
ment with the Hubbard model results as compared to the 
first approximation is noticeably improved as regards the 
skewing of the spectrum. 



c. Analytic approach The operator c\ a c 0: 



st 7. 



in Eq. |8| removes a spin singlet from the Heisenberg wave 
function and an ed 'singlet' is created: 



u 



n 



2(\de. 



\ed...)). 



The wave function of the final state, in the large-C/ limit 
can be written in the product formD: 



If) =R(H)®\x)® \<p)), 



(10) 



where \ip) describes L — 2 spinless free fermions on an L 
site lattice with twisted boundary conditioned (the holes 
represent the e and d), \x) is the squeezed wave func- 
tion of the remaining L — 2 spins with momentum Q = 
2wJ/(L-2), (J = 0,...,L-3) and \<p) = {\ed) + \de))/V2. 
The operator R = J^j e™ njrnjl corrects for the opposite 
sign of hopping of the d compared to e. Using this wave 
function, for Af(oj, k) we get 



AA(fc,c)=i^l6sin 2 ^ J2 



FnB, 



01 



xS(ui - E f + E GS )S k! p f . 



(11) 



The spin part gives F Q =J2\(Xq 2 |^1Xgs)I 2 > where the 
operator T removes a singlet from the first two sites of 
the spin wave function and the sum is over all states 
with momentum Q. This has been studied in detail in 
Ref. |l3], where it is found that «97% of the total weight, 
J2q Fq = (So ■ Si + j) (— > In 2 in the thermodynamic 
limit), is concentrated at Q = 0. This remarkable feature 
makes the calculation simple. 

The matrix element coming from the charge part of 
the wave function is B i — \(ipf\bobi\ipGs)\ 2 , where bi 



2 



annihilates a spinless fermion. First, we discuss the case 
V = 0. Then \ip) can be characterized by the quan- 
tum numbers of the two holes, l\ and l\ , with momenta 
Lk\ 2 — ^-^12 + Q> an< ^ the momentum and energy are 



Pf = -fcf -k% + Q=Y (-1? - ij + J) 



Ef = E GS + 2t (cos fcj 1 + cos &2 



(12) 



where we have neglected the small energy difference be- 
tween the exchange energies of the ground state and the 
final states of the order of t 2 /U, furthermore we assume 
that the ground state has momentum (i.e. the system 
has 2,6,10 etc. sites). Then B i reads 



Bqx = £2 sin " 



h,h _ h,h 
ft, j ft 2 



(13) 



Let us next consider the case with finite V: The effect 
of the nearest neighbor repulsion appears as attraction 
between the empty and double occupied sites , and leads 
to an extra — V(l — n 3 -)(l— Tij+i) term in the Hamiltonian 
describing the spinless fermions. Therefore we expect 
bound states (excitons) to appear in the spectrum. In 
the large U limit the factorized wave function remains 
an eigenstate even introducing the V terrrfl, supposing 
that V -C U. The two particle problem can be solved 
using the standard Bethe-ansatz 



\h,h) 



(14) 



where the holes are located at sites jx and ji. Since the 
system is translationally invariant, we can separate the 
momentum of the holes K (so that the total momentum 
is P f = -K+Q): 



cM(.L-h+ji) 



(15) 



From the twisted boundary condition aoj = —e l< ^a^L 
(we assume that the number of fermions and holes is 



±1, and the quantization 



even), v = - e i(KL/2-Q) 
LK = 2-kIk + 2Q follow, where Ik is an integer. Then 
the energy is simply (we introduce for compactness E' — 
E - U and J = u - U): 



K 

E = At cos — cosq, 
2 

and the secular equation for q becomes 



(16) 



„ Vcosq + 2tcos(K/2) 

2nl = Lq + 2 arctan — LJ -, (17) 

Vising 

where for v = 1 the quantum number / is integer, 
/ = 0,1,..., L/2 — 1, and half odd integer (/ = 
|, |, ...,(£ — 3)/2) for v = —1. Depending on the ratio 
V/[2t cos(K/2)], we can have either all the q's real, or, 
for sufficiently large ratio a bound state can appear with 



complex q. In the thermodynamic limit, L — > oo, the 
bound state is formed for V > 2t cos -|- with energy 



E' 



At 2 



K 



V -- C0S 2- 



(18) 



B i is now given by |a ,i| 2 / Y, 3l< j 2 \ a ji,h\ 2 and i aftcr 
a straightforward but tedious calculation, it reads 



B 



Wt 2 cos 2 %-E< 2 



oi — v 



L L(At 2 cos 2 f 



VE> + V 2 ) - VE> - 2V 2 



(19) 



At this point we have the weights to leading order in 
t/U. We have performed analogous calculations includ- 
ing also the term in the matrix elements, and the 
result obtained will be given later. Note that there are 
also contributions of the same order coming from the 
wave function correction, which would require going be- 
yond the factorized wave function. In order to investi- 
gate the relevance of these terms, we present in Fig. || a 
comparison of different approximations. Apart from an 
energy shift 0(t 2 /U) the dashed line is seen to be quite 
reliable, and is in fact better as far as the distribution of 
weights is concerned than is the dotted line, which shows 
that the operator corrections are more significant than 
those to the wave function. 

In the continuum limit we introduce g(K,uj') = 
L^2,j B x5(uj' — Ej), and we get 



9{K,J) 



m 2 cos 2 f 



,'2 



2tt At 2 cos 2 f 



At 2 



V 2 



Voj' 



- lO 

2 ,y ]s^'-E' exc ), 



6(4£cosf - 



(20) 



J g(k,Lu)du> — 1 and N{k,u>) now reads: 



M{k,u)ttyiQ S m 2 ^F Q g{k,u-U)+M inc {k,u), (21) 

where A/"inc(&,o>) contains the contributions from Q ^= 
and its weight is small (3%) compared to the main 
features given by the Q = part. 

The function g{k,uj) is shown in Fig. ^|. As previously 
alluded to, the spectrum consists of both continuum and 
sharp excitonic features. For V ~ t the exciton emerges 
from the continuum in the middle of the zone, and ac- 
counts for almost all of the spectral weight close to the 
zone boundary. For small momentum where the exciton 
is "inside" the continuum its' presence may still be seen 
as a strong enhancement of the intensity of the spectrum 
near the lower edge. Of course for V > 2t the exciton is 
present for all momenta. 

The effect of including in the operators can be 
given conveniently as 



g(k,u) g(k,uj) 1 - 



Auj + 2V 
U 



(22) 
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for the most important case of Q = 0. For V = 0, 
the t U~ 3 corrections only redistribute the weight and 
the sum rule fAf(k,oj)du = 161n2(t/J7) 2 sin 2 k/2 + 
0(t 4 /U 4 ) is satisfied. 

To conclude, we have determined the density response 
of a simple model for a correlated quasi-one dimensional 
insulator. The most surprising result is that excitons 
may appear near the zone boundary even if they are not 
present in the optical absorption (the k — > limit), de- 
pending on the parameters chosen. It would be interest- 
ing if these features can be observed experimentally. 

We acknowledge useful discussions with J. Fink, F. 
Gcbhardt, T. Giamarchi and P. Horsch. 
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FIG. 1. Af(k, ui) of a half-filled twelve site Hubbard model 
(full line) for V = t, U/t — 12, effective model with (dotted 
line) and without (dashed line) t /U corections in the density 
operators. The S functions are plotted as Lorentzians of width 
O.lt 
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FIG. 3. g(k,ui') for V = t and fc = (dotted), tt/2 
(dashed-dotted), 57r/6 (dashed) and it (solid line). The 
weight of the excitonic peak is not in the same scale as the con- 
tinuum. The inset shows the momentum dependence of the 
boundary the continuum and the exciton dispersion (heavy 
line) 
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FIG. 2. N(k,iu) for an 18 site cluster with U = 12t and 
V = It. The solid line is the "best" strong-coupling limit re- 
sult, where the exact wave functions and energies from the di- 
agonalization of H e fs as well as the n*- 2 - 1 correction is included, 
which are however omitted for the dotted line. The dashed 
line includes but uses the leading order wave functions 
and energies for the final states. 
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